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ABSTRACT 


A general  partitioned  transient  analysis  procedure  is 
proposed,  which  is  amenable  to  a unified  stability 
analysis  technique.  The  procedure  embodies  two  exist- 
ing implicit-explicit  procedures  and  one  existing 
implicit- implicit  procedure.  A new  implicit-explicit 
procedure  is  discovered,  as  a special  case  of  the  gen- 
eral procedure,  that  allows  degree-by-degree  implicit 
or  explicit  selections  of  the  solution  vector  and  can 
be  implemented  within  the  framework  of  the  implicit 
integration  packages.  A new  element-by-element 
implicit- implicit  procedure  is  also  presented  which 
satisfies  program  modularity  requirements  and  enables 
the  use  of  single-field  implicit  integration  packages 
to  solve  coupled-field  problems. 


TABLE  OF  CONTENTS 


Page 

INTRODUCTION  1 

LINEAR  COUPLED-FIELD  EQUATIONS  OF  MOTION  4 

IMPLICIT  INTEGRATION,  GENERAL  PARTITIONING  AND  EXTRAPOLATION  5 

EXAMPLES  OF  PARTITIONING  7 

STABILITY  OF  GENERAL  PARTITIONED  PROCEDURE  11 

IMPLICIT-EXPLICIT  PROCEDURES  IN  STRUCTURAL  DYNAMICS  15 

Element-by-Element  I-E  Partition  16 

DOF-by-DOF  I-E  Partition  17 

IMPLICIT- IMPLICIT  PROCEDURES  18 

Staggered  I-I  Partition  18 

Element-by-Element  I-I  Partition  20 

DISCUSSION  AND  CONCLUSIONS  22 

ACKNOWLEDGMENTS  24 

REFERENCES 25 


APPENDIX  A - DERIVATION  OF  CHARACTERISTIC  POLYNOMIAL  EQUATION  . . . . A-1 


APPENDIX  B - PROOF  OF  THEOREM  2 B-1 

APPENDIX  C - DERIVATION  OF  THE  STABILITY  Z-POLYNOMIAL  EQUATION  FOR  I 

ELEMENT  I-I  PARTITION  (44) C-1  i 


li 


INTRODUCTION 


The  transient  response  analysis  of  many  engineering  systems  involves  the 
formulation  of  semi-discrete  coupled-field  equations  of  motion  which  are 
then  solved  by  direct  time  integration  methods.  Such  coupled-field  equa- 
tions, for  example,  arise  in  the  modeling  of  structure-medium  interaction 
problems,  thermal-fluid  flow  problems,  and  structure-structure  interaction 
problems.  The  direct  time  integration  of  these  equations  presents  compu- 
tational difficulties  that  are  not  normally  encountered  in  single-field 
problems  or  in  structure-structure  interaction  problems  where  mesh  sizes 
are  more  or  less  uniform.  There  are  two  major  sources  of  difficulties, 
one  related  to  computational  aspects  and  the  other  to  implementation. 

First,  the  isolated  fields  may  have  vastly  different  response  character- 
istics. This  suggests  that  for  computational  efficiency  different  inte- 
gration algorithms  and/or  different  time  step  sizes  for  each  solution 
vector  should  be  adopted.  In  fluid-structure  interaction  problems,  for 
instance,  the  characteristic  response  time  scale  of  the  fluid-field  solu- 
tion vector  is  typically  much  slower  than  that  of  the  structure-field 
vector. 

Second,  the  bulk  of  existing  engineering  analysis  software  has  been  devel- 
oped primarily  for  the  treatment  of  single-field  problems.  One  is  thus 
motivated  to  make  use  of  existing  single-field  analysis  software  to  solve 
the  coupled-field  equations  of  motions;  this  Imposes  modularity  constraints 
on  software  architecture.  If  modularity  requirements  are  met,  the  direct 
time  integration  of  the  entire  equations  can  be  performed  either  by  se- 
quential or  parallel  execution  of  single-field  analyzers.  Solution  pro- 
cedures geared  to  this  operation  mode  will  be  termed  partitioned  transient 
analysis  procedures  (or  simply  partitioned  procedures) . As  shown  in  the 
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paper,  there  Is  a multitude  of  possible  partitioned  transient  analysis  pro- 
cedures for  each  coupled-field  system.  Therefore,  the  task  of  selecting  a 
"best"  procedure  can  become  formidable.  This  Is  because  the  analysis  of 
the  stability  and  accuracy  of  partitioned  procedures  is  much  harder  than 
that  of  single-field  solution  procedures.  Furthermore,  the  assessment  of 
I j the  relative  merit  of  competing  procedures  may  be  difficult  to  settle  If 

> different  analysis  techniques  are  employed  for  each  procedure. 

t Partitioned  procedures  for  structural  dynamic  problems  were  first  proposed 

^ by  Belytschko  and  Mullen  [1-2].  They  described  an  Implicit-explicit  (I-E) 

I partition  through  which  meshes  that  exhibit  high  (low)  frequency  response 

' . characteristics  are  treated  by  implicit  (explicit)  Integration  formulas. 

Hughes  and  Liu  [3-4]  proposed  an  alternative  I-E  partition,  which  appar- 
ently leads  to  a simpler  implementation  in  finite-element  programs. 

Few  systematic  studies  of  partitioned  procedures  for  specific  coupled-field 
problems  have  appeared  in  the  literature.  An  analysis  of  integration  pro- 
cedures for  the  coupled  sound  and  heat  flow  problems  by  Morimoto  [5]  appears 
to  be  the  first  study  on  the  subject.  Recently,  Park,  Felippa  and  DeRuntz 
[6]  studied  Implicit- imp licit  partitioned  procedures  for  a class  of 
fluid-structure  interaction  equations i Their  objective  was  to  treat 
the  symmetric  parts  of  both  the  fluid  and  structural  equations  by  implicit 
Integration  formulas  and  coupling  terms  by  extrapolation  formulas.  The  two 
primary  solution  vectors  were  the  scattered  fluid  pressure  on  the  "wet" 
structural  surface  and  the  structural  displacements.  Hence,  the  solution  of 
these  coupled-field  equations  was  obtained  by  a sequential  execution  of  fluid 
and  structural  analyzers,  which  led  to  the  term  "staggered  solution  proced- 
ures" . Unconditional  stabilization  of  such  implicit-implicit  procedures 
was  made  possible  by  exploiting  the  availability  of  the  fluid  radiation 
damping  term  in  the  equations  of  motion,  which  is  the  result  of  treating 
the  infinite  fluid  domain  by  a doubly  as3rmptotlc  approximation  [7].  The 
idea  of  introducing  damping  for  stabilization  was  subsequently  pursued  by 
Belytschko,  Yen  and  Mullen  [8]  in  their  study  of  an  implicit-implicit  pro- 
cedure for  a two-degree-of- freedom  model  problem  which  describes  certain 
fluid-structure  Interaction  problems. 
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Common  to  the  aforementioned  studies  Is  the  use  of  a specific  partitioned 
procedure  and/or  restriction  to  special  classes  of  problems.  The  present 
study  endeavors  to  consider  a general  partitioned  procedure  which  Is  amen- 
able to  a unified  stability  analysis  technique.  Motivation  for  undertak- 
ing this  study  Is  twofold:  It  eliminates  case-by-case  stability  analysis 
of  each  particular  partitioned  procedure,  and  It  provides  a common  ground 
for  evaluating  competing  procedures. 

As  a means  to  achieve  the  preceding  objectives,  the  present  paper  utilizes 
the  following  three  basic  formualtlon  steps.  First,  the  entire  seml-dls- 
crete  equations  of  motion  are  treated  by  Implicit  Integration  formulas. 
Second,  the  resulting  difference  equations  are  partitioned.  Finally,  coup- 
ling terms  that  are  transferred  to  the  right-hand  side  In  the  difference 
equations  are  extrapolated.  In  other  words,  the  first  step  completes  the 
temporal  discretization.  The  second  defines  computer  Implementation  as 
regards  matrix  configurations  that  have  to  be  treated  In  each  single-field 
analyzer.  The  last  step  determines  which  parts  or  portion  of  the  entire 
solution  vector  should  be  treated  explicitly.  It  will  be  shown  that  the 
present  formulation  can  embody  existing  partitioned  procedures  plus  hith- 
erto untried  ones  that  may  offer  algorithmic  and/or  Implementation  advant- 
ages. 
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LINEAR  COUPLED-FIELD  EQUATIONS  OF  MOTION 


We  consider  problems  governed  by  the  general  matrix  equation 

Ali  + Bu  + Cu  = f (la) 

where  u is  the  complete  solution  state  N-vector 


X y 

and  u , u . ..  represent  single-field  solution  state  vectors;  a superposed 
dot  denotes  time  (t)  differentiation  and  the  superscript  T on  a matrix  or 
vector  denotes  transposition.  Matrices  A,  and  £ represent  generalized 
mass,  damping  and  stiffness  matrices,  respectively,  and  f is  the  general- 
ized force  vector. 

The  term  "generalized"  is  used  to  emphasize  that  no  specific  form  of  mat- 
rices A,  and  £ is  required  for  carrying  out  the  basic  three  formulation 
steps  outlined  in  the  Introduction.  For  example,  if  A is  diagonal,  £ is 
null,  and  £ is  positive  definite,  equation  (1)  represents  a lumped-mass 
and  stiffness  system  that  describes  undamped  oscillatory  problems.  If 
the  entries  in  A corresponding  to  u and  those  in  £ corresponding  to  u 
are  zero,  equation  (1)  represents  a coupled  hyperbolic  (u  ) and  parabolic 
(u^)  system  that  describes  fluid-structure  interaction  problems  via  a 
doubly  asymptotic  approximation  [7].  Other  coupled-field  problems  can 
be  defined  by  appropriate  selection  of  entries  in  matrices  A,  £,  and  £. 
Thus,  the  application  of  the  basic  three  formulation  steps  to  equation  (1) 
can  cover  fairly  general  partitioned  transient  analysis  procedures  applic- 
able to  a wide  class  of  coupled-field  problems. 
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IMPLICIT  INTEGRATION^ 


GENERAL  PARTITIONING  AND  EXTRAPOLATION 


Implicit  linear  multistep  formulas  can  be  expressed  in  a compact  form  (see 
Felippa  and  Park  (9]  for  details): 


u 

-n 


= 6 u + h 

- n -n 


(2) 


where  h“  represents  the  contribution  of  historical  terms  and  6 is  a formula- 

dependent  generalized  step  size.  It  is  noted  that  each  field  solution  state 

vector  in  equation  (lb)  can  be  integrated  by  different  integration  formulas. 

Therefore,  equation  (2)  is  to  be  viewed  as  a symbolic  representation  of  these 

T X y T 

different  integration  formulas,  viz.,  for  u = (u  , u , ' 


6 = (6 


. ) '^ , we  have 
u 


),  and  similarly  for  the  historical  vector  h . Substitution 


X’  y’ 

of  equation  (2)  into  equation  (1)  yields 


n 


E u = g 
~ -n  “n 


(3) 


where 


E = A+6B  + 6 C 


g = f + (A  + 6 B)  h'^  + 6 A h 

®n  -n  ~ ~ -n  ~ -n 


(4) 

(5) 


Now  we  wish  to  solve  the  difference  equation  (3)  for  the  entire  coupled- 
field  equations  by  making  use  of  single-field  analysis  software.  To  this 
end  we  restrict  our  considerations  to  two-field  problems  and  consider  the 
following  general  matrix  partitioning: 


where 


E 

= E + E 

(6) 

'vx  ~y 

E 

= A+6B  +6^C 

(7) 

~x 

E 

= 6 B + (S^  C 

(8) 

Equation  (3)  is  rearranged  using  Equation  (6)  in  the  form 


E u = g - E u^ 
~x  -n  “n  ~y  -n 


where  u'  is  a suitable  extrapolator  for  u . 
-n  -n 


(9) 


We  call  the  formulation  as  given  by  equations  (2)  through  (9)  an  implicit- 
explicit  (I-E)  procedure  if  contains  non-zero  diagonal  entries,  and  an 
implicit- implicit  (I-I)  procedure  if  E^  contains  only  of f-diagonal  entries. 


We  use  here  the  term  partitioning  in  preference  to  splitting  for  the  fol- 
lowing reason.  The  term  "matrix  splitting"  is  often  associated  with  pro- 
cedures in  which  one  wishes  to  reduce  the  bandwidth  of  the  semi-implicit  itera- 
tion matrix  or  to  treat  directional-sensitive  matrix  coefficients  separately 
in  alternating  direction  methods  [10-12],  Thus,  in  splitting  procedures 
one  solves  for  the  entire  solution  state  vector  simultaneously.  On  the 
other  hand,  in  partitioned  procedures  one  solves  first  for  part  of  the 
entire  solution  state  vector,  which  is  then  used  to  solve  a portion  of 
the  remaining  solution  vector.  This  distinction  can  be  clarified  by  the 
following  examples  which  arise  in  structural  dynamics. 


EXAMPLES  OF  PARTITIONING 


For  reasons  discussed  later,  matrix  A is  not  partitioned  and  the  same  parti- 
tioning is  applied  to  ^ and  £.  Hence,  we  will  deal  only  with  matrix  £ in 
the  following  examples. 

Let  us  represent  the  matrix  £ as  a 3x3  block  matrix: 


c 

~xx 

<£xb 

0 

~bx 

4b 

~by 

0 

C . 

C 

~yb 

~yy 

and  the  corresponding  solution  state  vector  u by 


T . X b y.T 

u = (u  , u , u-') 


where  super  or  subscripts  x and  y designate  the  two  domains,  and  b desig- 
nates the  boundary  as  shown  in  Figure  1.  We  now  illustrate  some  useful 
partitions. 

Example  1:  Belytschko  and  Mullen  Node-by-Node  Partitioning  [1-2].  This 
partitioning  is  achieved  from  equation  (10)  by  grouping  the  boundary  con- 
tributions to  either  one  of  the  domains  and  choosing  C as: 

~y 


0 0 


0 0 0 


0 c c 

~ ~yb  ~yy 


This  partition  was  used  in  an  implicit-ejcpllcit  procedure.  Partition  (12) 

enables  us  to  obtain  the  entire  solution  vector  u in  two  separate  solution 

stages.  First,  the  solution  vector  u^  is  obtained  explicitly  as  A is 

” X b T ^ 

assumed  diagonal.  Then  the  remaining  vector  (u  , u ) is  obtained  by  solv- 
ing the  implicit  algebraic  equations  once  u^  is  available. 

Example  2;  Hughes  and  Liu  Element-by-Element  Partitioning  [3-4]. 

If  we  decompose  the  boundary  matrix  C. , in  equation  (10)  into  two  element- 

~bD 

level  matrices  according  to  the  two  groups  of  boundary-adjacent  finite  ele- 
ments, the  following  partitioning  for  C results: 

0 0 0 

C = 0 C.  (13^ 


'^bb 

c, 

•^y 

C , 

C 

/^b 

'^y 

in 

which  C 

where  C. , = C^,  + C?^,  , in  which  is  the  contribution  of  the  element- 
~ob  ~bb  ~ob  ~bb 

level  matrices  that  belong  to  elements  in  the  spatial  domain  x,and  simi- 
larly for  . 

~bb 


Example  3:  A New  DOF-by-DOF  Partitioning. 

A partitioning  suitable  for  working  at  the  degree-of-freedom  (DOF)  level 
results  if  we  take 


0 0 0 

^ ^ 

0 0 0 


0 0 c 

^ ^ ^yy 


It  will  be  shown  later  on  in  the  paper  that  this  partition  gives 
the  same  stability  condition  as  the  node-by-node  partition  (12).  Never- 
theless, this  DOF-by-DOF  partition  (14)  offers  cr.e  computational  advantage 
over  the  node-by-node  partition:  it  can  be  Implemented  within  a single-field 


implicit  integration  package  as  the  implicit-portion  of  the  stiffness 
matrix  remains  symmetric.  On  the  other  hand,  the  node-by-node  partition 
(12)  requires  that  integration  be  carried  out  in  two  separate  inte- 
gration packages:  explicit  and  implicit,  because  the  implicit-portion 
stiffness  matrix  is  not  symmetric. 


STABILITY  OF  GENERAL  PARTITIONED  PROCEDURE 


The  stability  of  the  general  partitioned  equation  (9)  can  be  examined  by 
seeking  nontrivial  solutions  for  f^  = 0 in  the  form: 

u = X u , (15) 

-n  -n- 1 

where  A is  the  solution  amplification  factor.  The  integration  formula  (2) 
can  now  be  expressed  in  terms  of  its  associated  characteristic  polynomial 
(see  Appendix  A) : 


[p(X)  - 6 a(X)]  u - 0 

— n— m — 


and  similarly  for  the  extrapolator  (A. 2): 


u^P^  = e(X)  u 
-n  -n-m 


(16) 


(17) 


where  the  subscript  m denotes  the  number  of  step  intervals  used  in  the 
formula  (2).  Substituting  Equations  (16)  and  (17)  into  Equation  (9)  one 
obtains  the  characteristic  system: 


J(A)  u = 0 

-n-m 


(18) 


where  J(X)  = [p^  A+6paB  +6p(a-X™+e)B 


(19) 


in  which  p stands  for  p(A),  etc. 
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Remark : It  is  important  that  velocities  be  updated  by  the  formula 


u 


6 u + h 
-n  -n 


(20) 


and  not  by  the  differentiation  fotnnula 


u = (u  - h ) / 6 
— n -n  -n 


(21) 


in  order  to  avoid  numerical  instability  (see  Appendix  A) . 


For  algebraic  convenience  we  Introduce  the  transformation: 


1 = 


1 + z 
1 - z 


(22) 


which  maps  the  solution-bounding  unit  disk  | X | < 1 onto  the  entire  left- 
hand  z-plane.  The  characteristic  equation  in  the  z variable  is 


(23) 


For  stability  the  roots  of  the  characteristic  equation  (23)  must  satisfy: 


Re(z^)  < 0,  i = 1,, 


, 2mN 


(24) 


where  N is  the  size  of  equation  (1)  and  m is  the  number  of  steps  entering 
in  the  integration  formula  (2). 


The  major  thrust  of  the  present  paper  is  to  treat  various  partitioned 
procedures  and  their  stability  conditions  by  a general  partitioned  pro- 
cedure (9)  and  a unified  stability  analysis.  To  this  end  we  will  limit 
ourselves  to  the  use  of  the  trapezoidal  formula  and  one  term  extrapolator 
in  the  sequel,  viz: 


12 


p(^) 

a(X) 


X - 1 


X + 1 


6 = 0.5  h 


e(X)  = 1 

so  that  equations  (19)  and  (23),  respectively,  reduce  to 

J^(X)  = (X  - 1)^  A + 6(X^  - 1)  B + 6^X  + 1)^  C 

^ 2 (26) 

+ 2 6(X  - 1)  By  + 46  XCy 

det  [(A  - 6 By  - 6^  Cy)  + 5 B 3 + 6^  c]  = 0 (27) 

The  stability  of  the  general  partitioned  procedure  (9) , when  the  trapezoidal 
integration  formula  and  one  term  extrapolator  are  used,  can  now  be  examined 
via  the  characteristic  z-polynomial  equation  (27)  subject  to  the  constraints 
(24).  As  basic  means  for  evaluating  the  stability  condition  (24)  we  cite 
the  classical  theorems  of  Routh-Hurwitz  and  Lyapunov  (see  Chapter  XV  of 
Gantmacher  [13]).  For  the  case  in  which  all  matrices  appearing  in  equation 
(27)  are  symmetric,  it  is  convenient  to  use  a theorem  of  Bellman  [14]. 

2 

Theorem  1 (Bellman).  If(A-6B  -6  C),B  and  C are  non-negative  defin- 

I — . I ^ . I ' /wy  ^ ^ 

ite,  and  either  (A  - 6 ^ - 6 ^)  or  C is  positive  definite,  then  equation 

(27)  has  no  roots  with  positive  real  parts. 


If  a particular  partitioning  gives  rise  to  an  unsymmetrlc  system  in  equa- 
tion (27),  i.e. , the  node-by-node  partition  (12),  the  theorem  of  Bellman 
is  not  applicable.  For  this  case  one  may  invoke  Siljak's  algebraic  criteria 
[15]  to  examine  the  positive  realness  of  equation  (27). 


Alternatively,  a perturbation  theory  such  as  the  one  expounded  by 
Vldyasagar  [16]  can  be  employed  to  assess  the  asymptotic  stability  of  equa- 
tion (27)  for  6 (0,»).  Stability  for  intermediate  stepsize  ranges  can 

then  be  determined  by  a model  problem  as  successfully  exploited  in  [6]. 


If  we  restrict  ourselves  to  undamped  cases,  i.e.,  ^ ^ in  equation  (27), 

a simple  criterion  can  be  derived  for  the  positive  realness  of  equation  (27) 
which  involves  unsymmetric  matrices.  To  this  end  we  rewrite  equation  (27)  as 


det 


A 

~x 


A -&^C^ 


6^  C 


= 0 


in  which 


(28) 


yx 

~y 


Sym[C^] 


C 

~y 


yy 

~y 


(29) 


In  equation  (29) , we  have  = C and  = C for  partition  (12) , and 

YY  yy  y y^ 

= jO  for  partitions  that  render  I-I  procedures.  We  now  state  a theorem 
for  equation  (28)  (see  Appendix  B) . 


Theorem  2.  If  ^ and(^  - are  non-negative  definite  and  C is  positive 

definite,  then  equation  (28)  has  no  roots  with  positive  real  parts. 


Remark;  Theorem  2 is  applicable  to  damped  cases  (jB  ^ 0^)  if  ^ is  a Rayleigh- 
type  [17].  This  is  because  jB  can  be  made  parts  of  A and  £ by  a commutable 
matrix  transformation  [18].  The  ensuing  stepslze  restriction  for  unsymmetric 
ally  partitioned  I-E  procedures  can  be  obtained  by  replacing  appropriate 
transformed  matrices  in  Theorem  2. 
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IMPLICIT-EXPLICIT  PROCEDURES  IN  STRUCTURAL  DYNAMICS 


In  the  previous  section  we  have  presented  a technique  for  determining  the 
stability  of  the  general  partitioned  procedure  (9)  through  the  concept  of 
the  solution  amplification  factor.  We  now  examine  the  three  partitions 
defined  by  equations  (12)  through  (14)  as  applied  to  structural  dynamics 
problems  when  the  trapezoidal  formula  and  one-term  extrapolator  (25)  are 
used.  It  will  be  shown  that  the  three  resulting  procedures  have  algorith- 
mic characteristics  identical  to  those  which  would  result  from  the  com- 
bined use  of  the  central-difference  formula  for  the  explicit  partitions 
and  of  the  trapezoidal  formula  for  the  Implicit  partitions. 

In  the  sequel  we  replace  the  matrices  in  equation  (1)  by  the  more  famil- 
iar notation: 


(A.  B,  £)  = (M.  D.  K) 


(30) 


Node-by-Node  I-E  Partition 

The  general  characteristic  z-polynomial  equation  (27)  for  the  partition 
(12)  becomes  for  D = 0: 


det  [(M  - 6^  K ) z^  + 6^  K]  = 0 


(31) 


Application  of  Theorem  2 to  equation  (31)  in  conjunction  with  equations 
(12)  and  (29)  indicates  that  this  procedure  remains  stable  provided: 


h < 2 / (m_„) 

BM  max 


(32) 
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( 


1 

■I 


! 


where 


('“’bm) 


is  the  largest  eigenvalue  of : 


max 


0 0 


0 K , 

~ ~yyj 


u 


(33) 


If  the  trapezoidal  and  the  central  difference  formulas  are  used  to  integrate 
X b T V 

(u  , u ) and  u , respectively,  it  can  be  shown  that  the  resulting  character- 
istic matrix  is: 


(A  - 1)^  M + (A  + 1)^  (£-&,)+ 


(34) 


It  is  noted  that  equation  (26),  which  is  obtained  by  specializing  the  general 
characteristic  equation  (19)  for  the  trapezoidal  formula  with  one-term  extra- 
polator,  is  identical  to  equation  (34). 

Element-by-Element  I-E  Partition 

For  the  partition  (13),  Bellman's  Theorem  is  directly  applicable  and  provides 
the  following  stability  condition: 


h < 2 / (a)„, ) 

HL  max 


(35) 


in  which 


w 


is  the  largest  eigenvalue  of : 


max 


(D  M u = 
HL  ~ - 


0 0 0 

° Srob  hoy 

P ’^yb  ~yyj 


(36) 


The  above  stability  condition  (35)  is  Identical  to  that  obtained  by  Hughes 
and  Liu  [3]  by  an  energy  stability  criterion.  It  can  also  be  shown  that 
the  combined  use  of  the  Newmark- family  predictor-corrector  pairs  for  the 
partition  (13)  gives  the  same  algorithmic  characteristics  as  equation  (34). 
The  presence  of  damping  can  be  accounted  for  to  give  the  same  stability  con- 
dition as  that  derived  in  [3]. 
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I 


t 


I ' 

( ■ 
{ • 


DOF-by-DOF  I-E  Partition 


As  the  partition  (14)  renders  a symmetric  system,  we  can  invoke  Bellman's 
Theorem  to  determine  its  stability.  It  can  be  shown  that  the  resulting 
stability  condition  of  this  partition  is  identical  to  that  of  the  node-by- 
node partition  (32).  The  algorithmic  characteristic  equation  is  obtained 
from  equation  (34)  by  taking  in  the  form: 


K 

~y 


0 0 0 
0 0 0 


(37) 


Remark  1:  Examination  of  the  three  I-E  partitions  demonstrates  that  regard- 
less of  the  particular  partition  selected,  the  computer  implementation  of 
the  present  procedure,  equations  (3)  through  (9),  can  be  accomplished  within 
the  framework  of  implicit  integration  formulas.  Hence,  any  need  for  explicit 
or  predictor  formulas  to  integrate  separately  the  explicit  partitions  is 
eliminated. 

Remark  2:  An  important  feature  of  the  general  partitioned  procedure  (9)  is 
that  it  allows  the  combined  use  of  the  element-by-element  and  the  DOF-by-DOF 
partitions  simultaneously.  This  is  possible  because  the  two  partitions  pre- 
serve symmetry  of  the  matrix  in  the  left-hand  side  in  equation  (9)  and  a single 
extrapolated  vector  is  applicable  to  the  coupling  term  in  the  right-hand 
: side  of  equation  (9).  For  example,  the  element-by-element  partition  can  be 

I 

I used  for  partitioning  different  elements  (shells,  beams,  low-  and  high-order 

I isoparametric  elements  for  solid  continua)  or  different  media  (fluid,  struc- 

Iture,  Soil  subdomains).  At  the  same  time,  the  DOF-by-DOF  partition  can  be 

applied  within  each  element  or  subdomain  to  integrate  the  DOF's  exhibiting 
low  (high)  frequencies  explicitly  (implicitly). 

I 

i 
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IMPLICIT-IMPLICIT  PROCEDURES 


The  need  for  Implicit- implicit  (I-I)  partitioned  procedures  arises  when 
the  coupled-field  equations  are  to  be  integrated  implicitly  by  sequen- 
tially employing  single-field  integration  packages.  We  will  now  consider 
two  possible  I-l  partitioned  procedures  that  may  be  useful  for  structure- 
soil  Interaction  and  structure-fluid  interaction  problems. 

Staggered  I-I  Partition 


Let  us  consider  the  partitioned  matrix: 


K 


~y 


0 

0 

0 


0 

0 

~yb 


(38) 


For  this  I-I  partition  the  solution-advancing  cycle  proceeds  as  follows. 

y 

First,  the  vector  u [see  equation  (11)]  is  implicitly  obtained  by  extrapo- 

by  X b T 

lating  u . Then  use  u to  obtain  implicitly  the  vector  (u  , u ) . 

Notice  that  the  characteristic  z-polynomial  equation  (28)  is  applicable  to 
this  I-I  partition.  Applying  Theorem  2 to  this  partition  in  conjunction 
with  equation  (38)  indicates  that  the  preceding  I-I  procedure  is  uncon- 
ditionally stable. 


This  is  in  striking  contrast  to  the  result  of  Belytschko,  Yen  and  Mullen 
[8]  from  their  study  of  a model  two-DOF  system,  in  which  they  proved  that 
the  preceding  partitioned  procedure  becomes  "unconditionally  unstable". 

We  offer  the  following  explanation  for  this  discrepancy. 
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The  difference  equations  employed  in  [8]  are: 


(M  + 6 K ) u 
~ -n 


= 2 (M-6  K)  u , - 

~ ~ -n- 1 


(M  + 6 K)  - 
~ ~ -n-2 


- 6^  K „<'■) 
~y  -n 


The  above  difference  equation  (39)  can  be  written  equivalently  by  a 
pair  of  difference  equations  and  the  integrator  for  velocity  as: 


(M  + 6 K ) u 
~ ~x  -n 


= (M  - 6 K)  u , + 2 5 M u , 
~ ~ -n-1  ~ -n-1 


Mu  = Mu  ,-6K(u  ,+u) 

— _n  — -n-1  ~ -n-1  -n 

On  the  other  hand,  our  procedure  (9),  when  specialized  to  the  case  of  the 
trapezoidal  formula  and  one-term  extrapolation  (35),  becomes: 


(M  + 6 K ) u 
~ ~x  -n 


(M  - 6^  K)  u , + 2 6M  u , - 6^  K u^^^  (41) 
~ ~ -n-1  ~ -n-1  ~y  -n 


Mu  = Mu  ,-6K(u  i+u) 

~ -n-1  ~ -n-1  -n 

Comparing  equations  (40)  with  (41)  one  sees  that  the  latter  is  equivalent 
(P) 

to  the  former  if  u = 0 in  equation  (41).  The  characteristic  z-polynomial 
equation  for  equation  (40)  can  be  derived  as 


det  (m  - is  6^  z^  - 4 Ky  z + 6^  K = 0 


from  which  one  can  establish  that  the  procedure  is  "unconditionally  unstable" 
for  the  model  problem  considered  in  [8]. 


Remark : If  the  extrapolator  in  equation  (39)  is  changed  to 


u'P’  . 2„  , 

-n  -n-1 


yn-2 


(43) 


it  can  be  shown  that  the  resulting  characteristic  equation  is  identical 
to  that  of  the  present  procedure  (41).  This  is  a direct  consequence  of 
the  well-known  property  that  the  order  of  accuracy  of  difference  equa- 
tions for  second-order  systems  maintains  that  of  the  integration  formula 
when  discretized  by  one-derivative  formulas  and  it  decreases  by  one  from 
that  of  the  integration  formula  when  discretized  by  two-derivative  form- 
ulas. Therefore,  from  the  viewpoint  of  our  procedure  (41)  the  procedure 
of  [8]  corresponds  to  neglecting  the  coupling  term  when  computing  u^. 


Although  the  present  staggered  I-I  procedure  (41)  which  is  a special 
case  of  the  general  procedure  (9)  is  unconditionally  stable,  it 
does  not  meet  the  modularity  requirement.  The  reason  for  this  is  that  the 
entries  in  the  matrix  in  equation  (38)  cannot  be  decomposed  into  the 
two  subdomains  as  shown  in  Figure  1.  The  remedy  to  this  shortcoming  is  an 
element-by-element  I-I  partitioned  procedure. 


Element-by-Element  I-I  Partition 


If  the  partitioned  matrix  ^ in  equation  (9)  contains  non-zero  diagonal 
entries,  the  solution  vector  pertaining  to  ^ is  necessarily  obtained  by 
an  explicit  process  as  shown  in  the  section  on  Implicit-Explicit  Proced- 
ures. One  can  overcome  this  deficiency  by  Implicitly  solving  two  sets  of 
difference  equations  for  both  partitioned  subdomains  (see  Figure  1). 

This  is  accomplished  by  the  following  I-I  procedure: 


(M  + 6^  K ) u _ = 

~ 'rx  -nt 

(M  + 6^  K ) u 
~ ~y  -n 

where  K takes  the  partition 


g - 6 K u , 
2v.i  ~y  -n- 1 


§n  - « -Sx  “nt 


(13). 


(44) 
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For  the  particular  case: 


K , and  E 
~x  ~y 


6^  K , the  sta- 

~y 


A = M,  E = M + 6' 

bility  of  the  element-by-element  I-I  procedure  (44)  can  be  analyzed  and  the 
resulting  z-polynomial  stability  equation  is  (see  Appendix  C): 


det 


(4M  + 2 6^  [k  K + K M ^ K 1 ) z^ 

\ ~ ~ Oy  ~ ~xj  ! 

+ 2 (k  K - K m”^  K ) z + 4 6^  K 

\~x  ~ ~y  ~ ~oc/  ~ 


(45) 


= 0 


The  first  and  the  third  coefficiert  matrices  are  positive  definite.  The 
second  matrix  (K  M^K  - K M'^K)is  antisymmetric  and  therefore  has 


~x  ~ ~y  ~y  ~ ~x' 

purely  Imaginary  eigenvalues  only.  Hence,  the  element- by- element  I-I  pro- 


cedure (44)  is  stable  by  Theorem  1. 


It  is  noted  that  the  preceding  I-I  partitioned  procedure  (44)  preserves  the 
program  modularity  of  each  single-field  implicit  integration  package.  As 
the  preceding  I-I  procedure  is  new  we  will  describe  its  computational  as- 
pects in  some  detail. 


X b T 

First,  equation  (44a)  is  implicitly  solved  for  the  vector 

b y ^T  solve  equation  (44b) 


utilizing  the  extrapolated  vector  (u  , , u^  , ) 

, -n- 1 -n- 1 

D y T X 

for  the  vector  (u^,  uj^)  by  utilizing  the  temporary  solution  vector  (u^^.» 

u^  )^.  The  final  solution  vector  to  be  used  for  the  next  step  integration 
X b y T 

is  thus  • If'  the  next  step  equation  (44b)  is  solved  first. 


The  integration  process  is  continued  by  solving  equations  (44a)  and  (44b) 
alternatively. 


From  the  foregoing  examination  of  the  computational  sequence  of  the  element- 


by-element  I-I  procedure  (44),  it  is  evident  that  two  implicit  solutions  are 


required  for  the  boundary  vector  u per  step.  Therefore,  the  price  paid  for 


realizing  the  computational  modularity  by  the  element-by-element  I-I  procedure 
is  one  additional  implicit  solution  of  the  boundary  vector  u^  as  compared 


to  that  of  the  staggered  I-I  procedure  (41). 


1 
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DISCUSSION  AND  CONCLUSIONS 


In  this  paper  the  complete  set  of  coupled-field  equations  of  motion  is  dis- 
cretized by  implicit  Integration  formulas.  A general  partition  is  then 
introduced  to  the  resulting  implicit  difference  equations.  Finally,  the 
coupling  terms  that  end  up  in  the  right-hand  side  as  unknowns  are  extra- 
polated to  advance  the  solution  cycle.  The  ensuing  stability  analysis 
for  the  general  partitioned  procedure  (9)  is  carried  out  via  the  concept  of 
solution  amplification  of  the  difference  equations.  We  now  summarize  the 
main  results  of  the  present  study. 

It  is  demonstrated  that  a unified  stability  analysis  is  applicable  to  all 
existing  and  some  new  partitioned  procedures  even  if  the  matrix  equations 
that  describe  the  stability  characteristics  of  the  partitioned  procedures 
are  not  simultaneously  diagonalizable.  This  applicability  is  in  fact  due 
to  the  equivalence  betv'een  the  well-known  energy  technique  [19]  and  the 
positive  definiteness  of  quadratic  polynomials  in  the  matrix  theory  that 
underlies  the  theorems  of  Routh-Hurwitz  and  Lyapunov. 

Our  partitioned  procedure  (9)  can  be  specialized  to  a general  implicit- 
explicit  procedure  that  embodies  those  of  Belythshko  and  Mullen  [1-2]  and 
Hughes  and  Liu  [3-4].  In  particular,  it  led  to  the  discovery  of  a new 
DOF-by-DOF  implicit-explicit  procedure  that  preserves  the  symmetry  of  the 
partitioned  matrices  and  allows  degree-by-degree  implicit  or  explicit  sel- 
ections of  the  solution  vector.  Furthermore,  our  implicit-explicit  pro- 
cedure can  accommodate,  within  the  framework  of  implicit  integration 
implementation,  the  concurrent  use  of  the  element-by-element  and  the 
DOF-by-DOF  I-E  procedures. 


t ; 


Our  procedure  (9)  becomes  the  staggered  I-I  procedure  (41)  when  the  parti- 
tioned matrix  does  not  contain  non-zero  diagonal  entries.  The  staggered 
I-I  procedure  is  unconditionally  stable  if  the  extrapolator  is  judiciously 
selected.  On  the  other  hand,  the  cause  of  unconditional  instability  of  the 
I-I  procedure  of  Belytschko,  Yen  and  Mullen  [8]  can  be  traced  to  an  improper 
choice  of  the  extrapolation  vector. 

A new  element-by-element  I-I  procedure  (44)  is  presented,  which  preserves 
the  program  modularity  of  each  single-field  integration  package  when  used  to 
solve  coupled-field  problems.  The  procedure  is  stable  and  its  only  computa- 
tional overhead  as  compared  to  that  of  the  staggered  I-I  procedure  is  the 
requirement  of  an  additional  implicit  solution  of  the  boundary  state  vector. 

So  far  we  have  proposed  and  analyzed  procedures  that  integrate  field  solu- 
tion vectors  with  the  same  time  step  size.  The  stability  analysis  of  pro- 
cedures which  integrate  field  solution  vectors  with  different  time  step  sizes 
is  much  more  involved  and  the  resulting  computational  gain  margins  are  in 
general  harder  to  assess.  Study  of  such  fractional-step  partitioned  proced- 
ures is  now  being  actively  pursued  and  will  be  reported  in  a separate  commun- 
ication. 
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APPENDIX  A 

DERIVATION  OF  CHARACTERISTIC  POLYNOMIAL  EQUATION 


The  linear  multistep  formula  (2)  can  be  expressed  In  the  form 


m 


(A.l) 


and  a stable  form  of  extrapolators  Is  found  to  be  [6]; 


“j  Vj 


The  historical  vector  h In  equation  (2)  Is 

-n 


(A.  2) 


(A.  3) 


provided  = 1,  which  we  assume. 


Equations  (A. 1-3)  can  be  expressed  In  terms  of  the  characteristic  value  A of 

the  difference  Equation  (10)  and  u , viz.: 

— n-m 


m “n  in  “ ® 

-n— m — n— m 


u^^^  * e(A)  u 

-n  -n-m 


h;;  - 6[a(A)  - a”]  u^,^+  [A™-  p(A)]  u„_^ 


(A.4) 


where  p(X) 


m 

5 


(A.4 

cont'd) 


The  partitioned  equation  (9)  can  now  be  expressed  as: 


a”  u^_^  - (A  + 6B)  Ka”  - p)  - 6(a“  - a)  u^_^] 


- « A(x’"  - P)  i - 6 (X“  - a)  (B  u + C u ) 
~ -n-m  ~ -n-m  ~ -n-m 


(A. 5) 


+ e E u =0 
~y  -n-m 


where  p,  a,  and  e stand  for  p(A),  a(A),  and  e(A),  respectively. 


The  auxiliary  equation  needed  to  eliminate  u from  (A. 5)  Is 

-n-m 


Apu  = 6Aau 
~ -n-m  ~ -n-n 


(A.  6) 


Substituting  the  equation  of  motion  (1)  Into  (A. 6),  one  obtains  for  f 


Apu  = - 5a  (Bu  +Cu  ) 

~ -n-m  ~ -n-m  ~ -n-m 


(A.  7) 


Combining  (A. 5)  and  (A. 6)  yields 

[p2  A + 6^  0^  + S^(a^  - a“  p + p e)  Cy 


I 


lO 


The  velocity  term  (A. 8)  can  only  be  eliminated  by  differentiation  formula 
resulting  from  (A. 4): 


u = (p/6o)  u 

-n-m  -n-m 


which  is  substituted  into  (A. 8)  to  yield: 


[p^  A + 6(7pBjj+  6^  0^  Cjj+  6 (pa  - x“p  + pe)  B, 


+ 6^(0^  - X®p  + pe)  CJ  u =0 

~y  -n-m 


(A.9) 


(A. 10) 


It  should  be  noted  that  if  the  differentiation  formula  (A.9)  is  used  in 
place  of  (A. 6),  the  following  equation  results: 


[p^  A + 6 pa  Bv  + 5(po  - x'*o'  + ae)  B, 


+ 6^  a^  C„  + 6^(a^  - X^a  + ae)  Cvl  u =0 
~x  ~y  -n-m 


(A.ll) 


This  has  roots  whose  magnitudes  are  greater  than  unity  in  general  for  con- 
sistent Integration  formulas;  thus  the  integration  process  is  unstable. 

The  reader  can  verify  this  assertion  for  the  trapezoidal  rule  and  one- 
term  extrapolator. 


APPENDIX  B 


PROOF  OF  THEOREM  2 


Let  us  consider  the  undamped  system  from  equation  (1) 


A o + C u * 0 


Equation  (B.l)  can  be  transformed  to  an  equivalent  form 


q + C q “ 0 


(B.l) 


where  q » L u 

^ mm 


(B.2) 


T * -1  -T 

A - L L and  C - L C L 
^ ^ 


The  characteristic  system  for  (B.2)  which  corresponds  to  equation  (28) 
can  be  derived  as 

z^Gq  + 6^Cq  » 0 


where 


r I 0 1 

_ 2-yx  i_,2cyy 
. ^ ^ ^ • 


(B.3) 


Equation  (B.3)  can  be  further  reduced  by  the  orthogonal  transformation 


I 0 


0 T 


(B.4) 


B-1 


so  that 


T C „ T = A = diag 
where  k is  the  size  of  C , . 

'^y 


H “k) 


Using  equations  (B.4)  and  (B.5),  one  transforms  equation  (B. 


O A A A 

z G V + 6 C V 


where 


G = 


C = 


_g2^T^yx  x-6^T^C^^T 


c 

~xx 


in  which  B • C T. 

~ ^y~ 


Case ; 


C and  C 


yx 


c . 

■^yx 


For  this  case  we  have 


« B^ 
= A 


Now,  equation  (B.6a)  can  be  expressed  in  the  form 


£1.  2 -1  s 

.h*  ^ & cuv 

lo  i-«y 


where  H is  a similarity  transformation  matrix 


H 


-1  T 
-A  B^ 


(B.5) 


3)  to 


(B.6) 


(B.7) 


(B.8) 


(B.9) 


I 


1 


9 2 2 

Note  that  If  ~ ^ ^ and  jC  are  positive  definite,  all  z^,  1 * 1,...,  N 

are  negative  real  as  the  similarity  transformation  preserves  the  eigenvalues 

of  original  system  (6.6a). 

Case:  = 0,  cf*  “ £.  • In  this  case,  we  have 

2 "I  2 2 

z G V + 6 C V = 0 


where 


1 


As  G is  not  dlagonallzable,  we  perturb  it  to  the  form 


(B.IO) 


.a 

I+6^el 


(B.ll) 


2 

where  0 < 6 e <<  1.  Substituting  the  perturbed  matrix  (B.ll)  into  (B.lOa) 
and  following  similar  steps  as  used  in  deriving  equation  (B.6a),  one  can 
obtain 


where 


^ -S  A _P  2 P A 

2 + a C r V = 0 

a I+S  el 


(B.12) 


(B.13) 


We  note  that  the  eigenvalues  of  a matrix  depend  continuously  on  the  elements 
of  the  matrix  (see,  e.g.,  Vidyasagar  [16]).  Therefore,  as  e 0,  equation 

p 

(B.12)  contains  the  eigenvalues  of  the  system  (B.IO).  Furthermore,  H is  a 

similarity  transformation  matrix;  hence,  if  £ is  positive  definite,  all 
2 ~ 

eigenvalues  z^,  1 = 1,...,  N in  equation  (B.IO)  are  negative  real.  This 
completes  the  proof  of  Theorem  2. 


B-3 


APPENDIX  C 


DERIVATION  OF  THE  STABILITY  Z-POLYNOMIAL 
EQUATION  FOR  ELEMENT  I-I  PARTITION  (44) 


The;  difference  equation  (44)  for  the  trapezoidal  formula  with  one-term  extra- 
polation becomes: 


(M  + 6^K  ) = (M  - 6^K)  u , + 26M  u , - 6^K  u„  . 

~ ~x  -n  ~ ~ -n-1  ~ -n-1  ~y  -n-1 

(M  + 6^K  ) u = (M  - 6^)  u , + 26M  u , - 6^K  u^°^ 

v-  ~yy-'  _n  ~ -n-1  ~ -n-1  ~x  n 


(C.l) 


from  which  one  obtains 


Mu^0> 

-n 


(M  + 6^K  ) u - 6^K  u , 

~y-'  _n  y -n-l 


(C.2) 


for  the  time  step  t = we  alternate  the  Integration  as 


(M  + 6^K  ) u^2  = (M  - 6^K)  u + 26  u - 6^K  u^ 

~ -n+1  ~ ~ -n  -n  '-x  -n 

9 22  ^0) 

(M  + 6 K ) u = (M  = 6 K)  u + 26  u - 6 K u''": 

~x  -n+1  ~ ~ -n  -n  -n+1 


(C.3) 


which  gives 


M u^?J  = (M  + 6^K  ) u . , - 6^K  u 
~ -n+1  ^x  -n+1  ~x  -n 


(C.4) 


Now,  substitute  equations  (C.2)  and  (C.4)  into  equations  (C.lb)  and  (C.3b), 
respectively,  and  subtract  the  two  resulting  equations  to  yield: 


M(u  - u ) + 6^K  u . , - 6^K  u = (M  - 6^K) (u^  - u 1) 
~'-n+l  -n  '-x  -n+1  ~y  -n  ~ ~ -n  -n-1 

+26  M(u  - u ,)  + 6^  m"^K  u - 6\  M"^K  u , 

~ -n  -n-1  ~ -n  ^ ^ ^ -n-1 


(C.5) 


The  velocity  terms  in  equation  (C.5)  can  be  expressed  by  the  integration  for- 
mula and  the  equations  of  motion  as: 


^(^n  - “n-l)  “ - «K(lin  Vl)  <^-6) 

Eliminating  the  term  (u^  - in  equation  (C.5)  by  equation  (C.6)  one  obtains 

+ S^i)  + «\(h„  + 2u„.,  + 


(C.7) 


+JSy!l'isJ  a„  U„.j)  . 0 

Finally,  transformation  of  X into  z via  equation  (22)  gives  the  desired  result, 
viz: 


det  [(am  + 2K  m"^K  + 2K  m”^K  )z^  + 26^(k  m‘^K  - K m‘^K  )s 


+ 46^K 


(C.8) 


It  is  noted  that  the  integration  sequence  which  starts  with  the  extrapolated 
-n-1  equation  (C.l)  changes  the  sign  of  the  second  matrix  in  equa- 
tion (C.8).  However,  as  the  eigenvalues  of  K m”^K  are  the  same  as  those  of 

_ '~X~ 

iSv~  following  equality  holds: 

y « 


T -1  T -1 

X K M K X - x^K  M K X 
- ~y~  ~X-  — ~x~  ~y- 


0 


(C.9) 


for  an  arbitrary  vector  x. 


Therefore,  the  I-I  partition  (44)  is  stable. 


